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SUMMARY 


A new approach to the finite element method wnich utilizes families of 
conforming finite elements based on complete polynomials is presented here. 
Finite element approximations based on this method converge with respect to 
progressively reduced element sizes as well as with respect to progressively 
Increasing orders of approximation. Numerical results of static and dynamic 
applications of plates are presented to demonstrate the efficiency of the 
method. Comparisons are made with plate elements in NASTRAN and the high- 
precision plate element developed by Cowper and his co-workers. Some 
considerations are given to implementation of the constraint method into 
general purpose computer programs such as NASTRAN. 


INTRODUCTION 


With the availability of general purpose computer programs, such as 
NASTRAN, at reasonable cost, utilization of the finite element approximations 
is common practice. In the conventional finite element method, a continuous 
structure is Idealized by discrete structural elements which are joined 
together at nodes. Structural characteristics are expressed in terms of 
nodal variables. Improvement of accuracy is generally made with respect to 
progressively reduced element sizes. If certain conditions are met, then the 
finite element approximation will converge to the true solution when the 
element sizes are reduced (Ref. 1), 

Unfortunately, reliable and practical error estimation techniques are t\ot 
yet available. In Important analytical computations it is usually necessary 
to complete two or more calculations of the same problem in order to establish 
the validity of t.e finite element model Itself. The calculations usually 
employ progressively refined finite element nets, altuough upper and lower 
bound estimates have been proposed also (Ref. 2). This is a tedious and 
costly process involving a considerable amount of duplicated effort. 
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Much interest has been shown in the development of high-precision finite 
elements so that better accuracy could be obtained with fewer elements* 
Incorporation of such elements into the NASTRAN program was reported to be 
under development (Refs. 3 and 4). 

Computational experiments as well as theoretical considerations have shown 
that, in terms of the number of variables needed to carry an analysis to a 
specified level of precision, the high-order or high-precision finite elements 
are more efficient than the low-order ones. This is particularly true in 
vibration and buckling analyses where eigenvalue problems must be solved in 
terms of the problem variables. An additional fact in favor of using fewer 
high-precision finite elements is that the number of necessary man-computer 
interface operations and the volume of data processing services are roughly 
proportional to the number of finite elements employed. 

In view of these findings it is logical* to explore potential benefits to 
be gained from the convergence process based on progressively increasing 
orders of polynomial approximation. In this convergence process, the finite 
element net is held constant and the order of polynomial approximating 
functions is varied. Existing error bounds such as that proposed by Fried 
(Ref. 5) indicate that the convergence rate will be exponential in this case, 
whereas the convergence rate is geometrical when the finite element sizes are 
reduced. It is noted that this error bound is valid only when the exact 
solution is sufficiently smooth and free from singularities. 

While there are many competing formats for stating finite element 
approximation problems, it was found that it is convenient to state the 
general problem as a quadratic programming problem. In this formulation, 
which will be referred to as the constraint method in the following 
discussions, the functional to be approximated (usually the potential energy 
expression) is written as a quadratic expression of the unknown coefficients 
of the polynomial approximating functions. The inter element continuity 
conditions and principal boundary conditions are stated as linear equality 
constraints. An advantage of this formulation is that all matrices that are 
necessary to define the numerical problem can be generated automatically for 
arbitrary orders of approximation. The finite elements so constructed will 
exhibit convergence with respect to reduced element sizes as well as with 
respect to increasing orders of polynomial approximation. Because of the 
latter type of convergence, it is unnecessary to reconstruct the finite 
element model when higher accuracy is sought. Additional advantages of this 
formulation are: (a) Since the unknowns are the scalar coefficients of the 

approximating polynomial sequences, it is not necessary to transform the 
variables and stiffness matrices from one coordinate system to another. 

(b) All finite elements can be made mutually compatible by specifying the 
appropriate connectivity through the constraint equations. This is a very 
Important feature of the new method because it permits consideration of 
structural stiffeners with greater ease than the standard finite element 
methods, (c) The new method will yield the axact solution when the exact 
solution is a polynomial with a degree less than or equal to the degree of 
the approximating polynomials, regardless of the number or orientation of the 
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finite element employed, (d) The accuracy of. solution and the computational 
efficiency are not sensitive to input numbering schemes (therefore, the 
method provides flexibility to the users in generating the structural models) . 
(e) The number of elements depends solely on the geometrical configuration of 
the structure to be analyzed, not on the desired degree of precision as in 
conventional analysis. Thus only the minimum number of elements sufficient 
for idealizing the structure needed to be defined. Other factors, such as 
existence of point loads and/or discrete supports, do not preclude the use of 
large elements. 

In the following, a solution technique utilizing the essential features of 
this formulation is discussed and applications of the constraint method are 
illustrated with numerical examples for static and dynamic analysis of plates. 
Comparison is made with results obtained by plate elements in the NASTRAN 
program and the 18 degrees-of-freedom plate element presented by Cowper et al. 
(Ref. 6). 


THE CONSTRAINT METHOD 


In the constraint method, the finite element approximation is treated as 
a direct energy minimization problem in which the minimum potential energy is 
sought subject to certain linear constraints. As in the conventional finite 
element method, the structure is idealized by discrete elements whose dis- 
placement characteristics are approximated by the polynomial functions defined 
over the element domains. Usually, the unknown variables are the coefficients 
in the assumed polynomials, although other variable definitions may be used 
also. The total potential energy is minimized with respect to these unknown 
coefficients subject t j constraints which ensure satisfaction of both inter- 
element continuity and cinematic boundary conditions. 

Detailed formulation of the constraint method has been presented elsewhere 
(Refs. 7 to 10). It will be outlined here as follows: 

The total potential energy tt Is obtained by assembling the element 
potential energies expressed in terms of the coefficients of the 
approximating polynomials as 

n ■ - 1/2 laj (Si |a| - (Zj | a | (1) 


In equation (1), (aj is a row vector, containing the polynomial coefficients 
and |a| Is transpose of t.aj ; [S] Is a symmetric, positive matrix containing a 
set of submatrices along Its diagonal; (ZJ Is a row vector associated with 
applied loading. This equation is treated as a quadratic objective function 
which is to be minimised subject to the following linear constraints: 
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( 2 ) 


IP] ja( . )r1 

where [P] and |r[ define the interelement continuity and the external boundary 
conditions. For homogeneous boundary conditions |r} is null. 

Several different algorithms can be used for solving the problem 
represented by equations (1) and (2), Most of these require separation of the 
independent variables from the dependent ones in the constraint equations. 

Then the problem can be reduced to solving a system ot simultaneous linear 
algebraic equations as explained in Appendix A, 

For structures subject to dynamic loading, inertia properties must be 
introduced in addition to the structural stiffnesses. In the case of free 
vibration, the equation of motion for the Kth element is expressed as 

K K * ^ 

where is the consistent mass matrix and ^ is the second derivative of 

|a| ^ with respect to time. The unconstrained equations of motion for the 
entire system are obtained by summation and can be written as 


[SI |a| + [M] |a| - 0 (4) 

After separating the independent and the dependent variables In the 
constraint set (see Appendix A), the unknown variables ja| can be expressed In 
terms of the free variables |a(.| , and the constrained equations of motion 
become 

[h'^ 1 [S] [H] ja^( + [M] [H] | a^ } - - [S] |h| (5) 

The matrix [H] and the vectors and |h | are defined In Appendix A. 

For homogeneous boundary conditions, |h| vanishes and equation (S) becomes 

(SI {aj + [M] -0 (6) 

where [S] and [M] are the constrained stiffness and mass matrices, 
respectively: 




ts] 

- 

[s] [H] 

(7) 

[M] 

' [h'^1 

[M] [H] 

(8) 


It Is noted that the eigenvalue problem associated with equation (6) Is 
relatively small since all dependent variables were eliminated. 


An Important feature of this formulation Is that all finite element 
approximation problems can be fully defined by the matrices [S], | Z | , [P], 

|r| , and [M] (for dynamic application) for arbitrary orders of approximation. 
These matrices can be generated automatically for any given problem. 


NUMERICAL EXAMPLES 


The efficiency of the constraint method Is Illustrated with examples for 
static and dynamic analyses of structural plates. A comparison Is made with 
results obtained using plate elements In the NASTRAN program and the hlgh- 
preclslon plate element presented by Cowper et al. (Ref. 6). Emphasis Is on 
the accuracy and convergence of the c^uatralnt method with respect to 
Increasing orders of approximation and using a minimum number of elements. 
Additional numerical results for static analysis of plates and shells can be 
found In references 7 to 10. 


Static Analysis 


The first example problem for static application is the simply-supported 
equilateral triangular plate (Fig. 1(a)) under uniform pressure q. The exact 
solution of this problem Is a 5th order polynomial (Ref. 11). 

The constraint method gave the exact solution when the 5th order 
polynomial was employed, and only one finite element was necessary. The 
results obtained by the 18 degrees-of-freedom hlgh-precislon element (also 
based on the 5th order polynomial) for various finite element layouts were 
presented in reference 6 for displacements and bending moments at the centroid 
of the plate. Figure 1(b) shows the layout given in reference 6 for N«1 and 
N«36, where N Is the total number of elements. Finite element layouts 
(N*25 and 100) used in the NASTRAN model are shown in figure 1(c). Due to 
symmetry, only one-half of the plate was considered. Results at the centroid 
of the plate are given in table 1. It Is seen tb .t 36 high-preclslon elements 
with 108 degrees-of-freedom (DOF) were needed to obtain precision to five 
significant digits whereas only 6 DOF were needed in the constraint method to 
achieve similar precision. The NASTRAN results were obtained by interpolation. 
Employing 100 CTRPLT elements with 166 DOF, 10 percent error was observed. 
Additional comparisons between NASTRAN plate elements and the constraint 
method are presented in figures 2 and 3 for the displaceskent and bending 
moment along the centerline of the plate, respectively. The NASTRAN 


555 



100-element model gave satisfactory answers for the displacements but only 
marginal accuracy for moment. Similar accuracy was observed for the bending 
moment alon'j the same line. 

The second example is a rectangular plate with two opposite edges simply 
supported, the third edge free, and the fourth edge fixed under uniform 
pressure q (Fig. 4(a)). This is an interesting problem because it comprises 
all common boundary conditions. Due to symmetry, only one-half of the plate 
was conslaered. Finite element layouts are shown in figures 4(b) and 4(c) for 
the constraint method and the NASTRAN model, respectively. The quadrilateral 
bending element CQDPLT was used in the NASTRAN model with 300 DOF. Result 3 
obtained by the constraint method were also reported elsewhere (Ref. 7). 

Rapid convergence was observed with respect to Increasing orders of approxi- 
mation. It was found that very good results were obtained for the 6th order 
approximation with 21 DOF (free variables). These are compared with the 
NASTRAN results in figures 5 and 6 for bending moments along a line in the 
middle of the rectangular pla:e. It is seen that correlation of the NASTRAN 
results for My with the exact solution is not as good as for M^. In this case 
the NASTRAN model overestimates the maximum My by about 50Z. It should be 
noted, however, that NASTRAN gave satisfactory results along the centerline 
of the plate. 


Dynamic Analysis 


The first example for dynamic application is a cantilevered triangular 
plate. Natural frequency of the plate was solved by the constraint method 
for various combinations of finite element layouts and orders of approximation. 
Results are given in table 2 together with the results obtained by the high- 
precision 18 degrees-of-freedom plate eletoent and the experimental data 
(Ref. 12). The results show that in the constraint method monotonlc 
convergence can be achieved by Increasing the orders of approximation as well 
as by reducing element sizes. It is noted that the DOF represent the total 
number of equations in the associated eigenvalue problems. Comparable results 
were obtained by the constraint method with fewer DOF. 

The next dynamic problem is the free vibration of a simply-supported 
square plate shown in figure 7(a). Two elements were used for one-half of the 
plate in the constraint method (Fig. 7(b)), and 200 elements in the NASTRAN 
model (Fig. 7(c)). Natural frequencies of the first three modes are presented 
in table 3. Monotonic convergence was obtained by increasing the orders of 
approximation in the constraint method. The NASTRAN results, presented in 
reference 13, are also given for comparison. It is significant that, the 
resulting number of DOF for the eigenvalue problem is much smaller in the 
constraint method. 
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IMPLEMENTATION 


Implementation of di^idld^intrilir following steps: 

algorithm given in Appendix A may be di>;iaeQ 

1. Define structural model 


b. 


d. 


f. 


(including order o£ epproxlnarlon Chat can 

-- inpur, 

Element and material propert - element ID and define 

Applied loads (referred ° '"^J^J^nftef or^oint ID if the 

point of application by reauired for distributed 

joint exists; only element ID is required 

load) Preferred to individual element 

f„1™deu« ircatll lor point aupporta; dellne elea»nt 
boundary number for line support) 


2. Generate and assemble matrix 


3. 


4. 


Unconatr.ln«l orJ.^Tf 

Onconatralnad jTpproxtaatlon) 

5” r i/i . i ^ <-trorr'«ri;priLtion, 

s:rrr ' conpatlglllty 

arJ «t.r ml boundary („„,1 „ eon.tant v.lu.) 

Enforced displacement vector jRJ mu 

Determine the rank of [P] acLmplished 

r,^S.‘;b'"rr tTri'" 

Constrained matrix generation 

a. <i..put. tt.»£or«tlo. «trlca. IH) («lu*tlo« (A7)).«d )h} 


a. 


c» 


d. 


e. 


b. 

c. 

d. 

e. 


(equation <^8)) y^tion (7)) 

ui'.'rr “l 
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5. Equation solver 




a. For static problem, solve [S] { * Izl - {r} for |a£.{ 

Compute ja^l (equation (A5)) and then |a| (equation (A4)). 
Separate |a| Into Kl . K«l, 2, . . , N, for each individual 
element. _ _ 

b. For dynamic problem, solve [S] | | + [M] { a^ f “ 0 

6. Output data processing 

a. Compute results for each element directly from the approxi- 
mating polynomials whose coef flci encs are determined in step 5. 

b. Users define the element ID and desired locations of output 
recovery: some default values may be provided. 


NASTRAN Implementation 


In executing these operations, step 1 requires some modification of 
NASTRAN procedures. In particular, the element compatibility data needed in 
constructing the constraint matrix [P] in step 2, and the options for 
specifying uniform line support conditions must be revised. Steps 2 through 4 
are new except that the cu)rent multiple-point constraints and enforced 
displacement in NASTRAN can be included in steps 2d and 2e. The current 
equation solvers in N.\STRAN may be used in step 5. New NASTRAN functional 
modules are also required for output data recovery in step 6, since the 
results are obtained directly from the approximating polynomials. 

It should be noted that finite elements generated by the constraint 
method can be combined with existing elements in NASTRAN if it is so desired. 
In this case, the unknown variables consist of both coefficients in the 
assumed polynomials and nodal variable components. The elements can be 
connected together by the constraint equations. 


CONaUDING REMARKS 


The constraint method is an efficient and cost effective approach to 
finite element approximations. It reduces modeling time significantly because 
fe%rer elements are needed. The structural model thus may ha generated faster 
and with fewer errors. The accuracy and computational efficiency are not 
sensitive to input numbering schemes, and remodeling is not required for 
greater accuracy. Results presented herein and those obtained in other test 
cases (Refs. 7 to 10) indicate that highly accurate results can be obtained 
by the constraint method at reduced computer costs. It is desiraole, however, 
to solve some larger problems to provide better comparisons between this 
approach and the current approaches to finite element structural analysis. 
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Efficiency of this approach may be further improved by the development of 
efficient algorithms for obtaining Such an effort is currently under- 
way at Washington University in St. Louis. 

Implementation of the constraint method into the existing general purpose 
computer program such as NASTRAN is considered f iasible and worthy of further 
investigation. 


% 

*5 
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APPENDIX A 


A SOLUTION ALGORITHM FOR THE CONSTRAINT METHOD 


The problem is to minimize the total potential energy (equation (Al)) 
subject to a set of constraints (equation (A2)): 

(Al) 
(A2) 


Min. TT = — laj [S] {a} - IZJ |a| 
Subject to: [P] |a| = {r[ 


We begin by selecting m linearly independent columns from [P] and renaming them 
[B]. Then equation (A2) can be written as . 


(B] |a|,| + iC] {ad = |r1 


(A3) 


where vector { a], | contains the variables associated with the linearly independent 
columns In [BJ , and {a^j contains the remaining variables in {af. Vector |a| is 
related to {ab[ and |a(.| by the following equation 




[T] 




(A4) 


In which [T] is the appropriate permutation matrix. 
From equation (A3) , we can write 


{ab| - [b‘^1 {r| - [b‘^ 1 [C] jad 


- 1 , 


(A5) 


Substituting equation (A5) Into equation (A4) , vector |a| can be expressed in 
terms of {ag | as 


where 


(a) - [H] |act +lhl 


tH] - (T] 


{h| - m 




(Ab) 


(A7) 


(A8) 
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Substituting equation (A6) into equation (Al) , the total potential energy it can 
be written as 


If = \ la^j [hT] [S] IH] \a^\ + 

iaj,j [HT] [S] {h} + I Ihj [S] {h} 

- lacJ [hT] {z} - Ihj {Z} (A9) 

Minimizing tt with respect to the elements of have 

[h'*^] [S] [H] jac} + [h’^ 3 ^[S]{h| - {z|^ =0 (AlO) 

Equation (AlO) represents a set of simultaneous algebraic equations. It is 
noted that the original n variables in |a^ were reduced to n~ra, where m is the 
rank of the constraint matrix [P]. When the boundary displacements vanish, 

|r[ is null. Equation (AlO) then becomes 

[H^l IS] [H] Jac} - [h’^I (z} = 0 (All) 

Once Ja^l is solved, Jat;} can be obtained from equation (A5) . 
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TABLE 1 SOLUTION FOR S.S. EQUILATERAL TRIANGULAR PLATE 


METHOD 


CONSTRAINT 

METHOD 




COWPER 


NASTRAN 


NO. OF 
ELEMENTS 


0ISPLACEMEN1 
AT CENTRO I I 

BENDING 
MOMENT 
AT CENTROII 

1 

6 

1.02880 

2.40740 

1 

3 

.617284 

1.08333 

36 

108 

1.02881 

2.40792 

25 

46 

.92 

1.78 

100 

166 

.99 

2.19 


TABLE 2 NATURAL FREQUENCY OF CANTILEVERED TRIANGULAR PLATE (STEEL, T 


FINITE 

ELEMENT 

UYOUTS 


CONSTRAINT METHOD 


T 

10" TYP. 
jL 


COWPER (REF. 6 ) 




6TH 

4TH 

5TH 

6TH 

15 

12 

20 

30 





36.6024 

36.5538 

36.5331 

36.5158 

139.187 

141.0743 

139. 3590 

138.9769 

194.499 

203.3356 
- -,.l 

194.0896 

L... ...I.,..- 

193.5854 


TABLE 3 NATURAL FREQUENCY OF THE SIMPLY-SUPPORTED SQUARE PLATE 


ORDER OF 


liWil 


MODE NO. 


CONSTRAINT METHOD 

NASTRAN 

4TH 

5TH 

6TH 

3RD 

6 - 

_ 12 

20 


.9298 

.9081 

.9069 

.9U56 

2.6972 

2.3962 

2.2782 

2.2634 

5.1170 

4.6325 

4.5474 

4.5329 


.9069 

2.2672 

4.5345 
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^ASTRAN-CTRPLT (1/2 PLATE) 

O 25 ELEMENTS 

A 100 ELEMENTS 

CONSTRAINT | 

I METHOD (EXACT) 


Y/L(X. 0) 

Displacement along centerline of the equilateral triangular plate 
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NASTRAN-CTRPLT (1/2 PLATE) 

0 25 ELEMENTS I 

1 

A 100 ELEMENTS 

— COHSTMINT 
1 METHOD (EXACT) 








Y/L(X.0) 


Figure 5*- Moment along centerline of the equilateral triangular plate 





























